#!/usr/bin/python
#-*- coding: iso-8859-15 -*-
import pyfits
import numpy as np

def append_field(rec, name, arr, dtype=None):
    arr = np.asarray(arr)
    if dtype is None:
        dtype = arr.dtype
    newdtype = np.dtype(rec.dtype.descr + [(name, dtype)])
    newrec = np.empty(rec.shape, dtype=newdtype)
    for field in rec.dtype.fields:
        newrec[field] = rec[field]
    newrec[name] = arr
    return newrec

class Catalogue():
	objs=[]
	hdulist=[]
	datatable=1
	ID_field='WDS'
	RA_DE_mode=1
	RA_field='RA'
	RAh_field='RAh'
	RAm_field='RAm'
	RAs_field='RAs'
	DE_field='DE'
	DEsign_field='DE-'
	DEd_field='DEd'
	DEm_field='DEm'
	DEs_field='DEs'
	VMAG_field='mag1'

	aux_datatable=0
	aux_ID_field=''
	

	def load(self,FITS):
		self.hdulist = pyfits.open(FITS)
		self.loadBasicValues()

		if (self.aux_datatable != 0 ):
			print "Loading aux values HDU:" +str(self.aux_datatable)
			self.loadAuxValues()

	def loadBasicValues(self):
		#print self.hdulist.info()
		#print self.hdulist[self.datatable].data.names
		#print self.hdulist[self.datatable].columns
		objects=[]
		fields=[('_ID','|S10'),('_RA','float32'),('_DEC','float32'),('_MAG','float32')]
		for field in self.hdulist[self.datatable].data.dtype.descr:
			if field[0]!='':
				fields.append(field)
		rectype=np.dtype(fields)
		print rectype
	
		for r in self.hdulist[self.datatable].data:
		        if (self.RA_DE_mode==0): 	
			#Valores de RA y DE directamente en radians
				print self.RA_field,self.DE_field
				ra=r[self.RA_field]*180/np.pi
				dec=r[self.DE_field]*180/np.pi
		        if (self.RA_DE_mode==1): 	
			#valores de RA y DE en campos separados hasta segundos
				ra=(r[self.RAh_field]+r[self.RAm_field]/60.+r[self.RAs_field]/3600.)*15
				if r[self.DEsign_field]=='+':
					signo=1
				else:
					signo=-1
				dec=signo*(r[self.DEd_field]+r[self.DEm_field]/60.+r[self.DEs_field]/3600.)
		        if (self.RA_DE_mode==2): 	
			#valores de RA y DE en campos separados hasta minutos
				ra=(r[self.RAh_field]+r[self.RAm_field]/60.)*15
				if r[self.DEsign_field]=='+':
					signo=1
				else:
					signo=-1
				dec=signo*(r[self.DEd_field]+r[self.DEm_field]/60.)

			record=[r[self.ID_field],ra*np.pi/180,dec*np.pi/180,r[self.VMAG_field]]
			record.extend(r)
			self.objs.append(np.rec.fromarrays(record,rectype))



	#Append to record aditional fields from auxliary table
	def loadAuxValues(self):
		aux_data=self.hdulist[self.aux_datatable].data
		IDs=map(lambda x:x[self.aux_ID_field],aux_data)
		objs=self.objs
		self.objs=[]
		fields=[('_ID','|S10'),('_RA','float32'),('_DEC','float32'),('_MAG','float32')]
		for field in self.hdulist[self.datatable].data.dtype.descr:
			if field[0]!='':
				fields.append(field)
		n_aux_fields=0
		for field in self.hdulist[self.aux_datatable].data.dtype.descr:
			if field[0]!='':
				f=('aux_'+field[0],field[1])
				fields.append(f)
				n_aux_fields=n_aux_fields+1
		rectype=np.dtype(fields)
		print rectype
		for s in objs:
			record=list(s.tolist())
			try:
				r=aux_data[IDs.index(s[self.ID_field])]
				record.extend(r)
			except:
				record.extend(np.zeros(n_aux_fields).tolist())
				pass
			#print record
			self.objs.append(np.rec.fromarrays(record,rectype))
			
		
	def filter(self,ra0,dec0,ra1,dec1,mag=20):
                s=filter(lambda x:(x['_RA']>ra0 and x['_RA']<=ra1) and  (x['_DEC']>dec0 and x['_DEC']<=dec1 and x['_MAG'] <= mag ),self.objs)
		return s
	

class WDSCatalogue(Catalogue):
    def __init__(self):
	print "loading WDS catalogue"
	self.datatable=1
	self.ID_field='WDS'
	self.RA_DE_mode=1
	self.RA_field='RA'
	self.RAh_field='RAh'
	self.RAm_field='RAm'
	self.RAs_field='RAs'
	self.DE_field='DE'
	self.DEsign_field='DE-'
	self.DEd_field='DEd'
	self.DEm_field='DEm'
	self.DEs_field='DEs'
	self.VMAG_field='mag1'
	self.load('WDS.fits')


class ngcCatalogue(Catalogue):
#Aladin Catalogue:
#catcat -fits  VII/118/ >ngc2000.fits
    def __init__(self):
	print "loading NGC catalogue"
	self.datatable=1
	self.ID_field='Name'
	self.RA_DE_mode=2
	self.RA_field='RA'
	self.RAh_field='RAh'
	self.RAm_field='RAm'
	self.RAs_field='RAs'
	self.DE_field='DE'
	self.DEsign_field='DE-'
	self.DEd_field='DEd'
	self.DEm_field='DEm'
	self.DEs_field='DEs'
	self.VMAG_field='mag'


	self.aux_datatable=2
	self.aux_ID_field='Name'	

	self.load('ngc2000.fits')

class HiparcosCatalogue(Catalogue):
    def __init__(self):
	print "loading Hipparcos catalogue"
	self.datatable=1
	self.ID_field='HIP'
	self.RA_DE_mode=0
	self.RA_field='RArad'
	self.RAh_field='RAh'
	self.RAm_field='RAm'
	self.RAs_field='RAs'
	self.DE_field='DErad'
	self.DEsign_field='DE-'
	self.DEd_field='DEd'
	self.DEm_field='DEm'
	self.DEs_field='DEs'
	self.VMAG_field='Hpmag'
	self.load('hip2.fits')

class crossCatalogue(Catalogue):
#Aladin Catalogue:
#catcat -fits  IV/27  >cross_index.fits
    def __init__(self):
	print "loading cross_index catalogue"
	self.datatable=1
	self.ID_field='HD'
	self.RA_DE_mode=1
	self.RA_field='RArad'
	self.RAh_field='RAh'
	self.RAm_field='RAm'
	self.RAs_field='RAs'
	self.DE_field='DErad'
	self.DEsign_field='DE-'
	self.DEd_field='DEd'
	self.DEm_field='DEm'
	self.DEs_field='DEs'
	self.VMAG_field='Vmag'

	self.aux_datatable=4
	self.aux_ID_field='HD'

	self.load('cross_index.fits')






def CostellationBounds():
#Data from stellarium 0.8 sources
		fi=open('./data/costellations_bound_20.dat')
		p=[]
		for l in fi.readlines():
			r=[]
			for f in l.split():
				try:
					r.append(float(f))
				except:
					r.append(f)
			pp=(r[0]*15,r[1],r[2])
			p.append(pp)
		fi.close()		
		return p

def CostellationFigures():
#Data from stellarium 0.8 sources
		fi=open('./data/constellationship.fab')
		p=[]
		for l in fi.readlines():
			r=[]
			for f in l.split():
				try:
					r.append(int(f))
				except:
				 	r.append(f)	
			p.append(r)
		fi.close()	
		return p



		

if __name__=='__main__':
#	C=WDSCatalogue()
#	ngc=ngcCatalogue()
#	C.readfile("./tycho/tyc2.dat.00")
#	C=HiparcosCatalogue()
	C=crossCatalogue()
	for s in C.objs:
		print s['_ID'],s['aux_Name']
#	print C.filter(0,0,45,45,plx=200)

#	CostellationFigures()
#	print C.stars[0][0],C.stars[0][1],C.stars[0][2],C.stars[0][3]
#	ngc=ngcCatalogue()
#	for t in ngc.ngc_names:
#	print t

#	print ngc.ngc
#	bsc=crossCatalogue()
#	print bsc.stars
